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ABSTRACT 


Title of Thesis: Investigation of a Model for Upward Flame Spread : 
Transient Ignitor and Burning Rate Effects 

Name of degree candidate: | Lee, Cheol Ho 

Degree and Year: Master of Science in Fire Protection Engineering, 1996 

Thesis directed by: Dr. James G. Quintiere, Professor, Department of Fire 


Protection Engineering 


Several studies have developed upward flame spread models which use 
somewhat different features. However, the models have not considered the transient 
effects of the ignitor and the burning rate. Thus, the objective of this study is to examine 
a generalized upward flame spread model which includes these effects. We shall 
compare the results with results from simpler models used in the past in order to examine 
the importance of the simplifying assumptions. We compare these results using PMMA, 
and we also include experimental results for comparison. The results of the comparison 
indicate that flame velocity depends on the thermal properties of a material, the specific 
model for flame lemgth and transient burning rate, as well as other variables including 
the heat flux by ignitor and flame itself. The results from the generalized upward flame 
spread model can provide a prediction of flame velocity, flame and pyrolysis height, 


burnout time and position, and rate of energy output as a function of time. 
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CHAPTER 1 INTRODUCTION 


Upward flame spread on vertical surface is a critical aspect of accidental fires 
because of its inherent high speed and potential consequences of fire growth to 
surroundings. Most of the principal researchers in the area of fire have devoted significant 
effort in trying to extend the knowledge on the mechanisms controlling flame spread and 


mass burning to represent this hazard and attempt to assess the relative contribution for a 


material. Here this research is interested in the effect of an ignitor, thermal inertia(kp c) of a 


material, and burnout during flame spread. 

Saito, Quintiere and Williams[1] developed a flame spread model which includes 
the relationship between flame height, pyrolysis height, and characteristic ignition time. In 
this model, flame height is controlled by heat released per unit mass of fuel consumed and 


mass loss rate per unit area, pyrolysis height depends on flame velocity, and characteristic 


ignition time is dominated by kpc of a material. They assume that the ignitor effect is zero, 


which means after ignition, mass loss rate is constant, that 1s steady burning. In other 
words, the ignitor effect, burnout effect, and unsteady burning are not included in the 
solution. 

The objective of this research is to develop transient flame spread model which 


utilizes the numerical solution based on the formulation outlined by Saito, Quintiere and 


Williams[1]. The model will be dependent on the different kpc values of a material. The 


model will be applied to a thermoplastic. Specifically, this research examines the model 
using polymethylmethacrylate(PMMA), as an example. 
The ultimate goal of the research is to examine the flame spread model, which 


includes the ignitor effect, burnout effect, and transient burning rate model performed by 


Hopkins[8], using the data obtained by some researchers[11,12,13] in the program and 
comparing the results with the experimental results of Orloff, de Ris, and Markstein[11]. 
The generalized results should provide more accurate predictions in terms of flame spread 
because it includes transient effects. Using the model we can predict the flame height, 
pyrolysis height, flame velocity, burnout position and time, total energy release rate at a 


specific time. 


CHAPTER 2 Derivation of The Exact Solutions of Flame 
Spread Model 


2.1 Review of “Upward Turbulent Flame Spread” by Saito, 


Quintiere, and Williams[1] 


2.1.1 Description of Spread Mechanisms 


Flame Spread occurs as a consequence of heating of the unignited portion of the 
fuel to a temperature at which vigorous pyrolysis begins. This heating is produced by 
convective and radiative heat transfer from the flames that bathe the fuel surface. Let x 


denote the vertical distance along the fuel surface, with x=0 at the base of the fuel, X=X, at 
the upper edge of the pyrolysis region and x=x, at the average height of the visible flame 


tip, as illustrated in Fig. 2.1. The heat transfer responsible for spread occurs in the region 


X2X_,. For steady-state burning at the base of a vertical wall, the energy flux q" to the wall 


has been found experimentally[2] to correlate with x/x, , and in a rough first approximation 


2 oe = ° 
for q" = qo = constant =~ 2.5 W/cm for O<x<x, and q"=0 otherwise, so that x isa 


good measure of the distance over which the principle heat transfer occurs. 
If this rough approximation is employed along with the further assumption that 


X_~ Xi) remains approximately constant during spread, then the upward spread velocity of 


pyrolysis front is 


2 
Vo= 44") (; — x.) / [mkpo(T, - T,)°1 (2.1.1) 


where k, p, c are the thermal conductivity, density and heat capacity, respectively, of the 
fuel, and T, and To are the ambient and ignition(or pyrolysis) temperatures of the fuel. 


Therefore, Equation (2.1.1) can be rewritten as 


Spi Schaant S | Ona 


a T, -- ir. 2 
where, cee q kpc Tie GG SOUS 


the characteristic ignition time tT for spread depends only on fuel properties, the ambient 


temperature and the level of the heat flux to the fuel from flame. As a simplification for 


describing time-dependent spread, we assume that Eq.(2.1.2) continues to apply with x, - 
X variable and that T remains an approximately constant time characteristic of upward 


spread. 


Fig.2.1 Illustration of the spread model 


2.1.2 Flame-Height Correlations 


Having hypothesized that the correlation of the heat-flux distribution with x/x, may 


lead to Eq.(2.1.2), we need an expression for X—- X, to obtain Mee By definition 


x0 = He + PV p(t at, (2ale3) 


where Xo is the value of x, at an initial time t=0, and 5 is the dummy variable of 
integration. Flame-height correlations are required for obtaining x, The total rate of 


energy release per unit length is the sum of 


Q' +qfrm"dx , (2.1.4) 
0 


where Q' is the energy release rate per unit length at the base of the wall, rn” is the rate of 


mass loss per unit area of the fuel, and q is the heat released per unit mass of fuel 


consumed. Flame-height correlations are of the form 


~ n 
x =k{ Origen dx), (2.1.5) 
f i, 0) 


where k, flame height coefficient, and n are constants. The flame height for wall flames is 
given such that k, = 0.067 (m?/kw’)!? and n=2/3, or approximately k, = 0.01 (m*/kw) 


and n=1[2],[{3],[4]. 


2.2 Exact Solution for n=1 

As a basis for describing upward spread we shall assume that the flame spreads 
after ignition(Q' =0) and the rate of mass loss per unit area(m" ) is constant. Following 
these assumptions and substituting n=1 into Eq.(2.1.5), the flame height can be rewritten 
as 


X= k, (q m" x) : (220) 


Substituting Eq.(2.2.1) into Eq.(2.1.2), Eq.(2.1.2) can be rewritten as 


Be (keqm" a0) Xy 


Ve |= 58S hee (2:27.23 
p dt t 


To derive x» Eq(2.2.2) can be rewritten 


ik (k qr" — 1) dt 
P _ ie ea (2.2.3) 
XD t 


Integrating Eq(2.2.3) x, can be obtained as 


(kqm" —1)t/T 
xy = ear ; (2.2.4) 


which means x, is increases with time(t). Therefore, substituting Eq.(2.2.4) into 


Eq.(2.2.2) the exact solution for n=1 is 


(kqm" —1)t/t 
dx (k qm" — 1)x ae 
Vic (2 DNR en eee (2.2.5) 
p dt 1 


2.3 Exact Solution for n=2/3 
Following the assumption Q' =0 and m" is constant and substituting n=2/3 and 
Bai(221-> )Anto.Eq(2.1-2); NE for n=2/3 can be expressed as 
ae Bia 
ah k (qm Xp) = Xy 


v= => = —— . (2.3.1) 
pa ordt - 


Unlike the case of exact solution for n=1, this case is required some steps to derive x, since 


Eq.(2.3.1) is non-linear. 
Let 


a. 1/3 (2i522)) 


1 os 
| oa ee (2.3.3) 
G 
and 
2 
ao = Ya fa fate (2.3.4) 
Substituting Eq(2.3.4) into Eq.(2.3.1), we have 
d¢ 2/3 
—_ = hn" ~ Brac. 


SATPA E 
n=k,(qm")" - ¢, (2.3.6) 


and differentiate both terms of this, then we have dc=— dn. Substituting these into 


Eq.(2.3.5), we have 


3 it Pao Ee 
T— = - eae 
tie (2.3.7) 
and 
d 
eal i (2.3.8) 
n oT 
After integration Eq.(2.3.8) we have 
eal SR (2.3.9) 


prt os 
No 


and substituting Eq.(2.3.2) into Eq.(2.3.6) and then substituting again Eq.(2.3.6) into 
FEq.(2.3.9) we have 


1/3 
Pp A - 1/3) 
1/3 


k-(qra")""” — x 
(2.3.10) 

2B 

k (qm ) — x 


and from Eq.(2.3.10) we have 


1/3 OAS: = 3 eee, 1/3 
Xp = k (qm ) —(e d{ k (qm ) Pay Lee (2.35.8) 


To show that X,, is some function of time, Eq.(2.3.11) can be rewritten as 


1/3 K2/3 1/3 Serta ace Sige 
X yy = k -(qm ) the x6 / kK (qm ete ke doe (273 2) 
and then »X, iS 
OU oe a Lo IS) et) St. 3 
xy = Ke (qm”") A Bs / k -(qm eae) ¢é oa (273113) 


which means x, increases with cubic time(t?). 


Therefore, letting he A in Eq.(2.3.13) and substituting A into Eq.(2.3.1), the exact 
solution for n=2/3 is 


dx k(qmn"A)~ ay 
a ee (2.3.14) 
Dias ci " 


CHAPTER 3 Derivation of the Flame Spread Model, 
and a Numerical Algorithms 


3.1 Integral Equation Formulation 


Since burning rate(m" ) is independent of the location of the element, the integral in 


Eq. (2.1.5) may be written as 


%p Po *» 
fm" dx = frm"dx +fm"dx , (3.105 
0 0 », 


where m" = m" (x(t), Liste t at X=X,(t). Since OSxSx,, Eq.(3.1.1) can be rewritten 


as 
*” *? 
a dxe= "rp (Xho? t) no + a: abe ae (3102) 


Eq.(3.1.2) shows that burning rate is related to the position of material, and all terms in 
Eq.(3.1.2) can be changed from the position to time since the position independent with 
time as shown Fig.3.1. Therefore, Eq.(3.1.2) can be rewritten as 


m0 t a 
eg dx = as Ce aie ah : (Si 


Substituting (dx /dt) oti Net into Eq.(3.1.3) and then substituting Eq.(3.1.3) into 


Egq.(3.1.2), Eq:(3:1-2) becomes 


10 


a t 
fm" dx = m"(t) x0 + ss (=t, Vt.) dt, (3.1.4) 


0 


Hence, substituting Eq.(3.1.4) into Eq.(2.1.5), the flame height(x,) become 


t 5 n 
Oued [Ope ms [m"t— ty V(t) dt, | (3.1.5) 


Xp = ky 


Therefore, substituting Eq.(2.1.3) and Eq.(3.1.5) into Eq.(2.2.2), the integral equation for 


flame spread is 


. 
Vp(t) =— r ' 


t Nn 
OQ +q [Ogg + rh"(t—t V(t) dt, i 


a Ee + JV Cp at | 


(3.1.6) 


xp(m) 


tp t time(s) 


Fig. 3.1 Illustration of pyrolysis front position response to the time 


11 


3.2 Numerical Solution 
To find the spread velocity with time, Vp(t), in Eq.(3.1.6), the integral equation in 
Eq.(3.1.6) should be solved. This study uses The Trapezoidal Rule[5] to solve the integral 


equation as a numerical method and an iteration process to find VO until convergence is 


satisfactory. 


3.2.1 Approximation Integrals by Trapezoidal Rule 


b 
The Trapezoidal Rule for approximating { f(x) dx is given by 


b 
b- 
IGS dx i [fx,) ni 2f(x ) SP PaEe 2i(x,_ 1) 35 Fx ,)] ’ (S220) 


To apply Eq.(3.1.6) to the Trapezoidal Rule let t=, 


I(j= J m'(t — t') Vp(t') dt’ , B22} 
1 i) 


and 


Lay J Vp(t') dt’ . (3233) 
a 0 


Following [n=1 > n=n+1] and [t’ =0 - t etm net atl Eq.(3.2.2) can be written as 


12 


Itt = f mt, — t)Vplt) at 
0 


+ 


: mtn - tv jt) = 
= a ar +m MS ay, — ty W p(tg) + 
4 may _ Lea JN: (t 
.. +m ae _ WY pn) ae 5) 
where h=t_.-t . 
n+1 n 
Defining @ as 
Mien: {1 luni ne 
One = t= —h 


13 


? 


(3.2.4) 


02 = te ote teen (3.2.5) 


where eo. = ty + (n)h, Eq.(3.2.4) can be rewritten as 
RR 
ihe oie 5 +m (G5 V (5) + 


"(04 DV n(ty..) 


w+ mG IV yn) + 5) (3.2.6) 
Therefore, Eq.(3.2.6) becomes 
I(t eee ¥ (a".V +l) aay ) (3:2 
Cea 2 ey ios Sep ja laos ce 
where m". = m"(@,) : o> aun - t. a es - (a-Dh , 
and Voi sas Vit) a ty + (a-I])h , t,=0 (ignition). 
Similarly Eq.(3.2.3) can be written following the process of above as 
‘hed 
E(t = RN pC ede 
2 n+l 0) 
an y (V..+V ) (3.2050 
2 Pi Pitl’ © = 


Therefore, following the Trapezoidal Rule, the integral equation, Eq.(3.1.6), can be written 


14 


as 


Vp(t) =—— tke /Q 4: q {th"(OX,, +1 © ah = [Xpo - LO} ; (3.2.9) 


3.2.2 The Solution of the Integral Equation by Iteration 
Assuming a new value(V ,(t)), which is in the Right Hand Side(RHS) in 


Eq.(3.2.9), to a previous value(V (t-1)), which is gotten from a previous step, we can get 
the new value(V (t)), which is in the Left Hand Side(LHS) and is not correct value. To 
find a real new value(V (0), some examination is needed like 


NAO mace) 


Error = v0 ate 8, (3.2014)) 


where € is a tolerance. If Error is greater than €, let V0 = yA] and then repeat the 


process until Error less than equal ¢. This will be shown later in computer program. When 


this condition is satisfied, we can get a new correct value(V ,(t)). 


15 


CHAPTER 4 Comparison of Exact Solution and Numerical 
Solution Using Computer Program 


A numerical solution is not exact since the solution comes from integral and 
difference approximations. We, however, will apply this numerical algorithm to a 
generalized flame spread model] that will be discussed later, therefore; we need to test its 
accuracy. To do this test, we shall compare the difference of results obtained from the 
exact and numerical solutions | 


The exact solutions used for testing are taken from Chapter 2 ; (1) x, « x, and (2) 
a sie . The exact solutions are given by Equations (2.2.4 )and(2.3.13). In both cases 


the mass burning rate per unit area is constant and the ignitor effect is zero. 


Since flame velocity depends on the differences between flame height(x,) and 
pyrolysis zone(x,,), we can predict the correlation as shown Figure 4.1. In both case, the 


velocity eventually becomes zero. In the case of n=2/3 the point, where flame height is 


equal to pyrolysis(x, = x): is earlier than the point in n=1. Therefore, the time in n=2/3, 
where flame velocity starts to decrease, is earlier than the time in n=1. This zero velocity 


point is an usual feature of both solutions, and may not be physical since t should decrease 


as the flame gets bigger. In any case they still form good tests for the numerical results. 
The numerical algorithm is programmed in Fortran. Also, the variables and data 
used for programs are is in Appendix A. The program is list in Appendix B. 
4.1 The Variables and Data used for Testing 
The material used for this comparison is PMMA and properties of this material are 


taken from Quintiere and Rhodes’s experiment[6]. The energy flux q" was already 


mentioned in Chapter 2. In steady state, the mass loss rate can be obtained from 


16 


ri" =———_—=_ , (4.1.1) 


where q" f is the heat flux from flame, and o is Stefan Boltzmann constant (5.67*10 


kw/m7k*), and L is heat of gasfication. The initial pyrolysis zone xe is selected as 0.3m 


and the input data are shown in Appendix A. All of input data used for the exact solutions 


(for n=1 and n=2/3) and the numerical solutions (for n=1 and n=2/3) are same. 


4.2 Programs for Testing 
Each program can be developed following the process described in Chapter 3 


respectively. In the program for the numerical solution, the tolerance(TOL) used for 
convergence 1s 10~ and the time step(H) is 1.0 second. These programs are shown in 


Appendix B. 


4.3 Comparison of Results for Testing 


Since flame velocity(V_) is related to flame height(x,) and pyrolysis zone(x_) as 
P f P 


shown Eq.(2.2.2), we just compare results of flame velocity obtained from the 
calculations. 

Figure 4.2 and 4.3 show the difference of flame velocity between the exact solution 
and the numerical solution at n=1 and n=2/3 respectively. The differences are negligible. 
Therefore, we can say the numerical solution procedure can be used in the generalized 


flame spread model with expected similar accuracy. 


uf 


xfo 


_ Xpo 
Xp 


Figure 4.1 The correlation between flame height and pyrolysis zone dependent on 


different powers (n). 
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Flame Spread Velocity, Vp (m/s) 
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Figure 4.2 Comparison of flame spread velocity for PMMA between exact and 


numerical solution for n=1 as a function of time. 
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Figure 4.3 Comparison of flame spread velocity for PMMA between exact and 


numerical solution for n=2/3 as a function of time. 
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CHAPTER 5 The Theory of Generalized Flame Spread 
Model 


We have discussed flame height and flame velocity after ignition and under constant 
mass loss rate. In general, flame height and velocity, however, can be affected by the 
ignitor, burnout, and transient burning rate. Therefore, we need a general model that 
includes the effect of an ignitor, burnout, and burning rate to analyze and predict a real fire 


situation. The model will be described below. 


5.1 Flame Height Calculations 


As shown Figure 5.1 flame spread can be separated with three parts. Figure 


xf(t) 


(A) (B) (C) 


Figure 5.1 Configuration of flame spread, (A)Before Ignition (B) After Ignition (C)After 


Bumout 


eal 


(5.1.A), (5.1.B), and (5.1.C) show the flame height effected by before ignition, after 
ignition and after burnout respectively. 
Flame height is solely due to the ignitor before an ignition occurs as shown Figure 
(5.1.A). Its flame height can be computed or experimentally determined according to its 
configuration[7]. For example, if its configuration is like a pool fire then 


= 0.23 oa = EO? D; 


g ; (53154) 


X fig 
where OF is the ignitor source(kw). If the ignitor is more like a line fire of width W 
against the wall, then 
See kpQi/W)" ) (5.1.2) 
where the correlation between k, and n is the same as before (Eq. 2.1.5). 
Figure (5.1.B) shows the flame height after wall ignites due to Qip and Q' up to 


burn out of the initial region ignited(t,, <t<t bX): At this situation flame height 
becomes 


x(0(= k-| Qi, / W) + Q' l (5.1.3) 


where Q' is the wall contribution. Figure (5.1.C) shows the flame spread after initially 


ignited burn out(t > t,(%_0))- At this time flame height can be written as 


X(t) = x, + k(Q)" (5.1.4) 


We can unify Figure (5.1) A, B, and C by introducing step functions 
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Nix lax = 0 
or 
TC)t 0 po 
Therefore, we can let 
igs) 


ers all (LG, Oe Qie : 


where, tx fig) is the burnout time at x=x fig’ also since Oe has a fixed duration time(At. ), 
ig 


the ignitor effect can be written as 


Qig = NGC )-D - (At, -1)- Qi, (5.1.6) 


Eq.(5.1.6) means that the ignitor can affect the flame height before burn out occurs at 


ES fig or before the duration time is achieved. 


5.2 Representation for the wall contribution (Q') and Burning Rate 


The wall contribution can be expressed as 


tp(x) tb(x) e 


Figure 5.2 Illustration of burning rate response to time 
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0 
Q = AH. - f m'(x) dx , 
0 


(Si) 


where AH. is the heat of combustion of a material (Before we used q in keeping with 


Reference[1]). 


As shown Figure 5.2 at position x the wall ignited or began to pyrolyze at time 


(x). This time corresponds to the time when x(t) =x. From previous work[8], we 


have an implicit formula for m"(t) at x, 


Gee 2K 


mh"(6)AH = Gas ors, _ “g lis 1): 
} A a) 2 2: AHL | By 3_ (8.8 
an Sata fe (X))) —In 
: GO MUT O8 6-5, 
where, On + : 
c(q" p ~ OT jp) 


a material constant for a specified flame heat flux, 


94 5(x) =  6axtt, (x) - te(x)) , depends on x 


AH _, is the heat of gasfication of a material, 


6 is thermal penetration depth of a material , 


t, is the time x (0) =x, t, is the time x(t) = EKA 
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(Giza) 


G.23) 


(5.2.4) 


(S255) 


Before Burnout 


LES Si 
| 
| 
| 
| 
| | 
| | 
CS eee on 
xfig xp 
(A) 
After burnout 
m" 
| 
| 
xb xp 
(B) 


Figure 5.3 Burning rate as a function of position 


ao 


The burning rate model assumes flame heating commences at t (x), and ignition occurs at 
(x). Each position x has its own burning history as shown figure 5.3. Note t,(x) is the 
time that x first experiences a heat flux due to the flame tip reaching x. The flame spread 


model assumes a uniform heat flux q" f from x to x, and zero heat flux beyond x é that is 


x>X,. Thus the flame spread model is 


dx Xp — X 
ee (5.2.6) 
Bp t At 
f 
Tt Tig = Ue 
where, At TE kpc = , aflame spread time, 
4 q f 


is constant for a given material. 


5.3 Representation for the wall contribution (Q') in terms of x 


As we discussed in Chapter 3.1, we need to consider the integral equation of 
Bq (2a baer 


x ti x, 


x (t) 
P 
I= — = f mx,t)dx= f m"dx+ f m"dx , (5.3.1) 
0 0 


x (t.) 
P 1g 


where, the meaning of the first term(I,) and second term(L,) of R.H.S is the burning rate 


in ignitor region and above ignitor respectively. Consider each term, 
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T= J m" dx? (5.3.2) 


since m" is constant over this region0<x < x(t.) = Xe Fq.(5.3.2) can be rewritten as 


Ly = m"-x (,,) =m" Xp - (5.3.3) 
And 
x= x(0)p, 
i= | ineate. (5.3.4) 
x= Mi) fig 


where, we need to convert to an integral over time. We recognize when t = t X= Xo and 


fi 


when t = to» X=X,. These relationships are the corresponding integral limits. 


Because m'"(x,t) = m"(@) , where 6 = bats (x), following the same process 
from Eq.(3.1.1) to Eq.(3.1.4), the relationship a = dx, /dt allows 


L, = jm - 1,00) “Vj (t,o) dt, (5.3.5) 


1g 


Since m"(x,t) is the burning rate per unit area at position x and at time t, if we know when 


x Started to pyrolyze 9), we can write down this value from our implicit formula, 


m'"(x,t) = m"(t — Ls m"(@(x)) , (5.3.6) 


Also, it is convenient to introduce tT where, 


oa 


t= t= tig : 56350) 
Therefore, L, becomes 
i= Vm (Oy Ve (teed 3. 
7 Jen, fF ons (5.3.8) 


8. AH | 5 = Stmy | Omas 
where, a) = = ____ In : 
6a | | Ps. 6-5, 
and r"(0) =|ay"- hy - a tu-TD| (Aree 


5.4 Burnout Effect 


We must limit m" due to burn out. For L, m" (t) is obtained from the formula m" ( 
6) where 6 =t-t (x__)=t-t, =t. Also burnout occurs after a duration 6 (x__) that is 
P eae ig b’ fig 
the duration for KEK the initial value. Hence the time for burnout is 


Vacate (5.4.1) 


tb es fig ig 


fig) = OL 


Or 


etl. (5.4.2) 


T(X fig ig 


fis? Ub fis’ > ‘bs 


As long as Tx ig) is greater than Tt this region continuous to burn. Hence we write 
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I, Sil oye NT, X Fig) — TT): Xfi ; (5.4.3) 


or 
Pen (0)s NO, Ago) —§(x)) - X fig > (5.4.4) 
where, 8 =T 
and O(x)= T- T(x). 
Similarly I, becomes 
Wy ifn@,co ~6) “1a"(8)| Vide (5.4.5) 


where, the burnout time is found by knowing the burning rate at that position X(T), 


G9 
m" = [ m"(6)d0, (5.4.6) 
0 


and m” is the burnable mass per unit area (g/m*). This can be found by the density(p ) of 


the wall fuel and its thickness(4 provided all the fuel vaporizes(m” = p 2). 


0 T (0) tT t 


Figure 5.4 The relationship between pyrolysis height and burnout position 


an 


The burnout position(x, (T)) can be found as 


X(T) = 0 for T < T, (0) : (5.4.7) 


which is before burnout of region OSXSX¢,, or 
X(T) = X(t’) for tT = T, (0) ; (5.4.8) 
which is after burnout of region x, , SxSx. 


Figure 5.4 shows the relationship between x, and X following Eq.(5.4.7) and (5.4.8) 


where Krone ‘ at T (0) and T( x(t’) is burnout time for x alte 
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CHAPTER 6 The Program and Results of Generalized 
Flame Spread Model 


Reviewing the theory of generalized flame spread model, a computer program can 
be developed. This program is more complicated than the program discussed in chapter 4 
since it includes the ignitor effect, transient burning rate, and burnout effect. From this 


program, we can obtain the pyrolysis zone(X,), the flame height(x,), the burnout 


position(x, ), the burnout time(T, ), the total energy release rate(Q), and flame velocity(V ) 


of a material at specific time(t). The program can be divided into four parts: (1) 


Declaration Part, (2) Initial condition, (3) Main Loop, and (4) Subroutines. Also the 
subroutines are separated into five parts, (1) ROOM to find steady penetration depth, (2) 
BURNOUT to find burnout time, (3) SEARCHB to find burnout position, (4) SPREAD to 
find pyrolysis zone, flame height, total energy release rate, and flame velocity, and (5) 
SEARCHF to find time of arrival of flame tip a material at specific time. 

As shown Figure 6.1 we can predict a typical result of generalized flame spread 
model, where time zero indicates the ignition time and x(1) are the initial values. The first 


drop and the second drop of x, occur when the ignitor is off and the burnout occurs 


respectively. The results of this model is in Appendix E 


6.1 Declaration part 
This part includes Input data which has the data of Material, Ignitor Characteristic, 


Heat Flux, Flame Height, and Computational Parameters, and includes the declaration of 


variables used for iteration, Computed Parameters,Maximum Number of Steps, and 
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0 


burnout 


alyep@lalieeha (epese 


Figure 6.1 The typical result of generalized flame spread model 
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Common Statement. Also this part includes the initialization part which has the data 
computed using input data of Material, Flame Height, and Computed Parameters. 

The material used for this example calculation is PMMA and properties of the 
material are picked up from Quintiere and Rhodes’ [6]. Energy output rate used for the 
ignitor characteristics, that is the size of the ignitor, and the flame height is picked from 
Back, et al.[9]. Also we chose the duration of ignition as 200s, that means the ignitor is 
turned off after 200s, and the width of wall heated as 0.5m. Heat flux from ignitor or wall 
flame is picked from the results of Williams, et al.[10]. All of data and variables used are 


shown Appendix C. 


6.2 Calculation Process 
This part is based on the theory of generalized flame spread model, and the 


program for this spread model will follow this process. 


6.2.1 Initial Conditions 
Hzet 
Tila 
c A) = te 
CALL BURNOUT (i, T(1), TG), m’”) - to find burnout time (t, (i)) at initial 


position 


x) = 0 


x) =x, +k, Qj, + Q). 
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where, Q atl (t,(1) - t()) - 11 (At, -t, -T(4)) Q 2 
and Q’, = AHe - rh" (i,i) , 
CALL ROO™ (i, t(i), i, t(i), 4,(4), t,@, m"(Li)) > to find burning rate rh" (i,i) 


x (1) = Xe 


X pW) — x) 
Vp(i) = 
At f 


Qi) = (Qing + QV) w 


6.2.2 Main Loop 


We have computed i values and we seek i+1 


T(itl)= T(i)+h 
Since we need T GQ) to begin and it depends on knowing X (+1) we must “guess” by using 


the previous value. This only affects the calculation of oe (i+1) and should not present a 


significant error as long as h is small. We will compute tT itl) after finding x (+1). 


a" (+1) = ih x) =x A) , the flame is heating x aig Sx) from time 0 


Or 


T.Ci+1) = TC) 
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CALL BURNOUT (i+1, t(i+1), T,G+1), T (i+]), m”) - to find burnout time 
(t, G+1)) at time i+] 


x (i+1) S10) TUt+1) < Tt 1) ,rfegionO<x< ne has not burned out 


Or 


x (+1) = x) ; T(i+1) 27, (1) 
CALL SEARCHB ( T(i+1), x, (i+1)) 
This finds x C+1) = X (4) where TG j= TO+1). 
CALL SPREAD (i+1, Vat), x (i+1), x (i+1), Q(i+1)) 


CALL SEARCHF (i+1, x G+), T,G+1)) 


6.2.3 Subroutine ROOTM (i, T(i), J, TG). 7 (). T @). mm" Gj) 
This finds burning rate m" (i,j). T(i) is the current time and LAC) is the time 


corresponding to xp(t)=x. From Eq.(5.2.2), Eq.(5.2.3), and Eq.(5.2.5), 


m'(i,J) = = Vi oe fey rei 


Vv 


LlGas () = (i)) | Ak | 


6-5) {7 - 10 
5 =8,— (5,—5,,) exp z Pre ee 


tT 


s 
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Where, q' = q's — oT. i , het flame heat flux, 


1g 
2kL : 
6. - p> aaa wee a material constant for a specified flame heat flux, 

e(q" - OT i.) 

2 

* 6. AH, 
and GS seer Sand burn time constant, 

6a 


and 
04 5(X) =  Ox(t, (x) ~ te(x)) ' 


To find 6 the program will use an iteration loop. For first guess to iterate, a previous value 


is used. 


6.2.4 Subroutine BURNOUT G, TG). 7 D. TS 1).m”’) 


This finds burnout time (T, ). T i) is the time of arrival of flame tip at position 


X=X, (i) and m” is a burnable mass per unit area. From Eq.(5.4.6), 


TH) 
m" = f m"(jdt - 
0 


Using Trapezoidal Rule to solve the integral, this equation can be rewritten as 
n-1 


m" = > [m"(ij) + m"G+1,j)] . 


i=] 
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The burnout time is when the integral >m” where the last i=n-1 and n=i+1. That is, t(n) = 


(n-1)h. Therefore, TG) = (n-1)h. Here, we need subroutine ROOTM to find burning 
rate m"(i,j). That is, 

ROOTM (i, T(i), j, T) TQ), eo, mi" (1,j)) 

ROOTM (i+1, t(i+1), j, LUE TG)» eo aeiny)_., 


where is the value of burnout time in subroutine ROOTM. Since we are integrating up to 


the burnout time, we can put this value(c°) as a big number. 


6.2.5 Subroutine SEARCHB (T(i), x, (@)) 


This seeks j such that TQ)= T(i). 

Dog =i 

IF [t,0) < T(i)] Then Continue Do Loop 
Else x,(1) = x9), for 7,9) 2 TG) 


Return 


End 


6.2.6 Subroutine SPREAD (j. V (j). xj). x (j). Q(j)) 


This find the velocity of the pyrolysis front, Vp(j), the flame tip position, x,(j), the 


pyrolysis front position, x0); and the total energy release rate, Q(j). 
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Qi =n (F,()-7H)- 0 t,t, 10) Q,, 


OT aH oe X pig * 
We need initial guess like ek, = VjG-1) to find a real V0) by an iteration loop. 


FIG) = 2/h [ FIQ-1) + m"G, 3-1) ViG-D + m'Gj) V@I . 


j-2 
where, FIG-1) => Sy (rh"G.K)V (k) dam j,k +1) ECAR 


is can be gotten following the same step described in Chapter 3.2. We need call ROOTM 


to find burning rate at that time and position. 


Q, =AH, FI(j) 
xD) = kp Qigt Q_ + Qn) +O - 


Following the same way discussed in Chapter 3.2 to solve integral for Xo» 
x) =C + (h/2) VG) . 


2 


h 
where, Gis X fig itt ti iE (V)@) + Vj (k+ 1)) + V 57 »| 


k= 1 


x-(j) — x_,G) 
Vi) eee 
, Ate 


Vp@ = VpG-D| 


Error = - e. 
V 
pY) 
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When this condition is satisfied, we can get a new correct value(Vp(j)). 
6.2.7 Subroutine SEARCHF (j, x (0. % (Wd 
This finds the time when the flame tip first reached X=X (j)- 


Do. k= 1.j 
IF [x,(k) < x0)] Then Continue Do Loop 


Else x(k) = x0) 
t,@) = (k-1) h 


Return 


End 


6.2.8 The Program of Generalized Spread Model 


The program follows the process described above. To make the program simple 
we use common statements, and the time step(H) is 1.0 second. 


Subroutine ROOTM is called by subroutine BURNOUT and SPREAD to find 


burning rate as described above. However, burnout time, Th has a different value when 


ROOTM is called by these subroutines. For example, once ROOTM is called from 


BURNOUT, we put Th with an “infinity’(big) value, however, Th is put with its true 


value, that is found in BURNOUT, when ROOTM is called from SPREAD. Therefore, 


the subroutine ROOTM is only used to find 5, and whenever the subroutine ROOTM is 


required to find burning rate, this program writes down the equation of burning rate after 
the Call ROOTM statement. 


The Program of Generalized Spread Model is shown in Appendix D. 
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CHAPTER 7 Comparison of Results 


The velocity of flame spread is related to the pyrolysis front position of material, 


x (i), the flame tip position, x {i and the characteristic ignition time,T, that is affected by k 


oc, as shown Eq.(2.2.1) and Eq.(2.2.2). In this section we compare the relationship 


between x and x, and the relationship between Me and x, of the exact solution for n=1, 


n=2/3, and the generalized flame spread model with the results that others found for 


PMMA. 


7.1 The Properties Used for Comparison 


Orloff, de Ris, and Markstein [11] reported, in their experimental study, k = 
0.64*10°° cal/em°C, p = 1.19 g/cm?, andc = 0.50 cal/g°C, respectively. Therefore kpc 


is 0.654 kW’s/m*C*. These values were assumed constant over the temperature range 


from ambient(20°C) to ignition (363°C) and under heat flux 25 kw/m*. They measured the 
burning rate of PMMA during upward flame spread finding it varied from 7.2 to 12.0 


g/em*s. Their initial condition was taken as 0.02 m. We chose an average burning rate of 


9.6 g/m*.s and Secon 0.02 m for our “exact constant burning rate solution” comparisons. 


The variables and data used for the comparison are in Table F.1. 


Mitler and Steckler[12] used the LIFT value derived for kpc of PMMA(1.02 


kW?s/m‘C”) in their study. We will only use this value for this comparison, and the other 


properties are same with Orloff, et al data. The variables and data used for the comparison 
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are in Table F.2. 


In case of the generalized flame spread model, we use the kp c from the LIFT data 


and 1kW of energy output rate(Qi,) which is calculated to make Xone m. The other 


properties are same with the data described in Appendix C. Using the different ignition 


temperature(T’,), also, we try to find the effect of ignition temperature(T) on flame 
spread. The ignition temperatures(T;,) used for this are 180°C and 363°C that come from 


J. Quintiere and B. Rhodes’ [6] and L. Orloff, et al.[11] respectively. 


Table 7.1 shows the different kpc values used for the comparison. 


Table 7.1 


The kpc properties of PMMA used for the comparison 


PROPERTIES 


Orloff et al 7.63*104 1190 0.654 


Generalized Flame 0.346*10° 1180 Ze) 1.02 
Spread Model 
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7.2 The Relationship between xX, and x, 


Orloff, de Ris, and Markstein measured the relationship, labelled experiment result 


between x, and x, finding a best fit of 


0.781 
= 1.95 Xp : GEFEN 


as Shown in Figure 7.1. Using their properties(q=25 kW/m’, m" =9.6 g/m’ .s) and 


Eq.(2.2.5) we can find the relationship, labelled exact solution as 


Xe = 24 xy . s10Gn= le (7.252) 


and 


Kp = 2.59 ta See prey: (7.2.3) 


Similarly we can also find the relationship using the PMMA LIFT data. The flame 
height relationships, however, are same with those of Orloff, de Ris, and Markstein 


because the same q and m" are used. The results of the flame velocity for these two data , 


however, will be different since they have the different value of kpc. These results will be 


shown later. 
Delichatsios, Mathews, and Delichatsios[13] found the relationship using 0.052 as 


a flame height coefficient,k,, 
(7.2.4) 
as shown Figure 7.2. 

The properties used for the relationships are in Table 7.2 
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Table 7.2 


The properties used for the relationship between 


pyrolysis height, x, and flame height, Xe 


PROPERTIES 


Source 


1 0.01 
Orloff et al 
213 0.067 
& a 


UNITS 
kW/m kW/m? g/m? S 
m 
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Xf vs. Xp 


Xf=1.95*(Xp**0.781) 


Xf(m) 


Flame Height, 


0.0 Ws 0.4 0.6 0.8 1.0 


Pyrolysis Height, Xp(m) 


Figure 7.1 The relationship between flame height and pyrolysis height for PMMA by 
Orloff, de Ris, and Markstein 
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Xf vs. Xp 


Aft=2, 00 (spr 705667) 


Xf(m) 


Flame Height, 


0.0 (6) 0.4 0.6 0.8 1.0 


Pyrolysis Height, Xp(m) 


Figure 7.2 The relationship between flame height and pyrolysis height for PMMA by 


Delichatsios, Mathews, and Delichatsios 
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7.3 The Relationship between Vp and x, 


Orloff, de Ris, andMarkstein measured the relationship called experiment result 


between Me and x, with 


V_= 0.00441 x Came k 
Pp Pp 


(Toa) 


as shown Figure 7.3. Using the relationship of Eq.(7.2.2) and Eq.(7.2.3) and substituting 


these equation into Eq.(2.2.2) we can find the relationship of yo and x, in the exact 


solution as 


V -7=_0.01448 x= Gator n=)". 
p p 


and 


0.667 _ 


V_ = 0.0103(2.59 x X_) , forn=2/3. 
P Pp p 


Similarly we can also find the relationship of Ms and x, for LIFT data as 


V.nt=10.009283: x=forn=t-- 
p p 


and 


0.667 _ 


V_ = 0.00663(2.59 x Xx.) 5. for n=2/3.. 
Pp p Pp 


(F332) 


(733) 


(7.3.4) 


(Tdas) 


The properties used for the relationships are in Table 7.3, where T is ignition time. 
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Vp vs. Xp 


0.008 


Vp=0.00441*(Xp**0.964) 


o 0.006 
= 
a 
> 
= 
= 0.004 
2 
a 
> 
® 
= 
© 0.002 
LL 


Pyrolysis Height, Xp(m) 


Figure 7.3 The relationship between flame velocity and pyrolysis height for PMMA by 
Orloff, de Ris, and Markstein 
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Table 7.3 


The properties used for the relationship between 


flame velocity, Vy and pyrolysis height, Xo 


PROPERTIES 


Source 


l 
Orloff et al 0.654 
ING 
1 
LIFT 1.02 150.799 
2/3 


Oe OS 
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7.4 The Programs used for Comparisons and Results 
The programs used for comparison are in Appendix G. The relationship of 


pyrolysis height, x? and flame height, x,, and flame velocity, var and pyrolysis height, 
X > are shown Figure 7.4 - 7.13. Figure 7.8 is the result of the comparison of flame height 


and pyrolysis height between the exact solutions and the experiment and the generalized 
flame spread model. These curves in figure 7.8 show the effect of the different flame 
height coefficient and power to the flame height. Figure 7.13 is the result of the 
comparison of flame velocity and pyrolysis height between the exact solutions and the 


experiment and the generalized flame spread model. These curves in figure 7.13 also show 


the effect of the different kpc and ignition temperature(T,, .) to the flame velocity. 
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Xf vs. Xp 


3 
= 2 
2: Xf=2.4*Xp 
= 
— 
@® 
<= 
E 
& 1 
LL 
0 
0.0 0.2 0.4 0.6 0.8 1.0 


Pyrolysis Height, Xp(m) 


FIGURE 7.4 The result of flame height vs. pyrolysis height used exact solution for 
n=1 with Orloff, de Ris, and Markstein data. 
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FIGURE 7.5 The result of flame height vs. pyrolysis height used exact solution for 
n=2/3 with Orloff, de Ris, and Markstein data. 
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Xf vs. Xp 
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FIGURE 7.6 The result of flame height vs. pyrolysis height used exact solution for 
n=1 with LIFT data. 
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FIGURE 7.7 The result of flame height vs. pyrolysis height used exact solution for 
n=2/3 with LIFT data. 
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Xf vs. Xp 
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(1) xf = 2.59 * (xp**0.667) , Exact solution for n=2/3, based on constant m" 
(2) xf = 2.4 * xp, Exact solution for n=1, based on constant m" 


(3) xf= 2.01 * (xp**0.667), Delichatsios, Mathews, and Delichatsios, based on 
constant m”" 


(4) xf= 1.95 * (xp**0.781), Experiment by Orloff, de Ris, and Markstein 
(5) Generalized flame spread model based on transient rh" 


FIGURE 7.8 The comparison of flame height vs. pyrolysis height. 
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Vp vs. Xp 
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FIGURE 7.9 The result of flame velocity vs.flame height used exact solution for 


n=1 with with Orloff, de Ris, and Markstein data. 
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Vp vs. Xp 
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FIGURE 7.10 The result of flame velocity vs.flame height used exact solution for 


n=2/3 with Orloff, de Ris, and Markstein data. 
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Vp vs. Xp 
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FIGURE 7.11 The result of flame velocity vs.flame height used exact solution for 


n=1 with LIFT data. 
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FIGURE 7.12 The result of flame velocity vs.flame height used exact solution for 


n=2/3 with LIFT data. 
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Vp vs. Xp 


0.02 


Vp(m/s) 


0.01 


Flame Velocity, 


Pyrolysis Height, Xp(m) 


(1) Vp = 0.0103*(2.59xp**0.667 - xp) , Exact solution for n=2/3 with Orloff, de Ris, 

and Markstein data, kpc = 0.654 kW?s/m4C2, rm" = const, T,,=363°C 

(2) Vp = 0.01448*xp, Exact solution for n=1 with Orloff, de Ris, and Markstein data, 
kpc = 0.654 kW2s/m4C2, m" = const , Tjg=363°C 

(3) Vp = 0.00663*(2.59xp**0.667 - xp) , Exact solution for n=2/3 with LIFT data, 
kpc = 1.02 kW2s/m4C2, m" = const, Tg=363°C 

(4) Vp = 0.009283*xp, Exact solution for n=1 with LIFT data, 
kpc = 1.02 kW2s/m4C2, m" = const , Tjg=363°C 

(5) Generalized flame spread model, kpc = 1.02 kW?s/m4C?, transient m" Ti =] 80°C 

(6) Vp = 0.00441*(xp**0.964), Experiment by with Orloff, de Ris, and Markstein 
7.2g/m2.s <m"<12g/m2.s,Tj=3630C 

(7) Generalized flame spread model, kpc = 1.02 kW?s/m4C2, transient m" 11} =363°C 


FIGURE 7.13. The comparison of flame velocity vs.flame height. 
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CHAPTER 8 The Effect of Thickness and the Ignitor 
on Flame Spread 


Using the generalized flame spread model with kp c=1.02 kW?s/m‘C? and the 


properties described by Quintiere and Rhodes [6] in Appendix C, we try to find the effect 
of thickness and the ignitor on flame spread in this section. A study on the effect of 


thickness and the ignitor include variations of thickness(mm): 0.1, 0.5, 1.0, 3.0 ; ignitor 


duration(s) : 30, 60, 120, 480 ; Q' ; oo : 10, 25, 50 or correspondingly Xo (mi) 


0.5, 1.0. Figure 8.1 - 8.3 show for the very thin material and low durations of the ignitor, 
the flame will never reach 5 m. But as these parameters are increased, propagation occurs 
and at faster speeds. Figure 8.4 shows the critical values of the parameters on propagation 


to 5 m. It is clear that all of these factors play a critical role in propagation. 
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FIGURE 8.1 Time to reach 5 mas a function of material thickness and ignitor duration 


at 10 kW/m for the ignitor. 
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Time to reach 5 m as a function of material thickness and ignitor duration 


at 25 kW/m for the ignitor. 
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FIGURE 8.3 Time to reach 5 m as a function of material thickness and ignitor duration 


at 50 kW/m for the ignitor. 
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FIGURE 8.4 Estimated critical values for propagation to 5 m. 
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CHAPTER 9 CONCLUSIONS 


Using the formulation outlined by Saito, Quintiere and Williams we developed a 
numerical algorithm that was checked with exact solutions for n=1 and n=2/3. 

We developed a general solution which included transient burning rate m" (t) and 
used a specific function for thermoplastics, and included burnout and ignitor effects. 


The comparisons illustrate the effect of model(that is, n=1, 2/3 and m" is constant) 


and value of kpc. Flame height and flame velocity in case n=2/3 have been found to be 
greater than those in case of n=] at initial, then these are switched. Also the effect of kpc 


has been found that the bigger kpc, the lower flame velocity. We also compared to 


experimental results, showing the general solution is in better agreement than the simpler 
analyses which assume m" is constant. That is because the case of transient burning rate 
requires less energy than that of constant burning rate. 

Future work should show a range of results from the general model to illustrate 
more clearly the role of the ignitor included duration time and heat flux and burnout related 


to thickness of material. It should be noted that any m" (t) function can be used in the 


general algorithm. 
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APPENDIX A VARIABLES AND DATA USED FOR TESTING 
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Properties Name of Data Units 
Variable in 
Program 


Ee 0.432103 | kWimk | 


(m5/kW2)1/3 


TABLE A.1 The Variables and Data used for Testing 
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APPENDIX B PROGRAMS FOR TESTING 
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REAL KAY, QFLX,QFLX0, XPO, CONST, Q, K, DEN, C, TP, TA 

REAL TAU, XP(0:1000), XF(0:1000), VP(0:1000), M(0:1000) 

DATA KAY, QFLX, QFLXO,XPO, CONST, Q, N / 0.01,25.0,25.0,0.3,5.4,0, 1/ 
DATA K, DEN, C, TP, TA, PI / 0.000432, 1190.0, 4.12, 375.0, 25.0,3.14/ 


OPEN(UNIT=11, FILE='CASE1.DAT', STATUS = 'UNKNOWN') 


TAU = ((PI*K*DEN*C)*(TP-TA)) / (4*QFLXO**2) 
DO 100 I=0,1000 
M(1) = CONST 
XP(I) = XPO*EXP((KAY*M(I)*QFLX-1)*/TAU) 
XF(I) = KAY*((Q+(QFLX*M(D*XP(1)))**N) 
VP(I) = (XF(1) -XP(1))/TAU 
WRITE(11,444) I, XP(1), XF(), VP) 


100 CONTINUE 

444 FORMAT(10X, I5, 5X, E11.6, 5X, E11.6, 5X, E11.6) 
STOP 
END 


Program 1 Exact Solution for n=1 
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REAL KAY, QFLX,QFLX0, XPO, CONST, Q, K, DEN, C, TP, TA 

REAL TAU, XP(0:1000), XF(0:1000), VP(0:1000), M(0:1000) 

DATA KAY,QFLX,QFLXO,XPO,CONST,Q,N /0.067,25.0,25.0,0.3,5.4,0.0,0.667/ 
DATA K, DEN, C, TP, TA, PI / 0.000432, 1190.0, 4.12, 375.0, 25.0, 3.14/ 


OPEN(UNIT=11, FILE='CASE2.DAT’, STATUS = 'UNKNOWN’) 
TAU = ((PI*K*DEN*C)*(TP-TA)) / (4*QFLXO**2) 
DO 100 I=0,1000 
M(1) = CONST 
A = EXP((-1*1/(3*TAU)) 
B = 1 - (XPO**(1./3.))/(KAY*((M(I)*QFLX)**(2,/3.))) 
CS ASB 
De (14 0)**3 
E = (KAY**3)*(M(I)*QFLX)**2 
XP(1) = D*E 
XF(I) = KAY*((Q+(QFLX*M(1)*XP(1)))**N) 
VP(1) = (XF(I) -XP())/TAU 
WRITE(11,444) I, XP), XF(I), VP(1) 


100 CONTINUE 

444 FORMAT(1X,15,1X,E11.6,1X,E11.6,1X,E11.6) 
STOP 
END 


Program 2 Exact Solution for n=2/3 


fi 


PARAMETER(NMAX=1000, NAR=1000) 

REAL KAY, QFLX,QFLX0, XPO, CONST, Q, K, DEN, C, TP, TA 

REAL Gl, G2, G3, N 

REAL H, VP(0:NAR), M(O:NAR), XP(NAR), X(NMAX), SUM1(0:NAR), SUM2(0:NAR) 
DATA KAY, QFLX, QFLXO,XPO, CONST, Q, N / 0.01,25.0,25.0,0.3,5.4,0, 1/ 

DATA K, DEN, C, TP, TA, PI / 0.000432, 1190.0, 4.12, 375.0, 25.0, 3.14/ 

DATA TOL, H /0.0001,1 / 


OPEN(UNIT=11, FILE='CASE3.DAT', STATUS = 'UNKNOWN') 


45 
200 
50 
201 


100 
444 


SUM1(0) = 0.0 
SUM2(0) = SUM1(0) 
TAU = ((PI*K*DEN*C)*(TP-TA)) / (4*QFLXO**2) 
VP(0) = (1/TAU)*((KAY*(Q+(QFLX*CONST*XPO))**N)-XPO) 
DO 100 I = 0,999 
M(I) = CONST 
M(I+1) = M() 
VP(I+1) = VP(1) 
ITERATION LOOP 


K=0 

X(1) = 0.0 

K = K+1 

IF(K .GT.NMAX) GOTO 45 
G1 = M() * VP) 
G2 = M(I+1) * VP(I+1) 
G3 = VP(I) + VP (I+1) 
SUM1(I+1) = SUM1(1+G1+G2 
SUM2(1+1) = SUM2(1)+G3 


VP(I+1) = (1/TAU)*(KAY*(Q+QFLX*(M(I+1)*XPO+ 
((H/2)*(SUM1(I+1)+G1+G2))))**(N)-(XPO+((H/2)*(SUM2(I+1)+G3)))) 
X(K) = VP(I+1) 


FIRST ITERATION. NO CONVERGENCE TEST 
IF(K .EQ.1) GOTO 10 
ERR = ABS((X(K) - X(K-1))/X(K)) 
IF (ERR .LE. TOL) GOTO 50 
GOTO 10 
WRITE(11,200)'DID NOT CONVERGE’, VP(I+1), ERR 
FORMAT(AI6, E8.2, 5X, E8.2,) 
WRITE(11,201)'CONVERGED AT’, K,'-TH ITERATIONS WITH ERR="", ERR 
FORMAT(A13, I5,A24,3X,E8.2) 
WRITE(11,444) I+1, VP(I+1) 
CONTINUE 
FORMAT(10X, I5, 5X, E11.6) 
STOP 
END 


Program 3 Numerical Solution for n=1 
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PARAMETER(NMAX=1000, NAR=1000) 

REAL KAY, OFLX,QFLXO0, XPO, CONST, Q, K, DEN, C, TP, TA 

REAL Giy GZ, G3, N 

REAL H, VP(0:NAR), M(O:NAR), XP(NAR), X(NMAX), SUM1(0:NAR), SUM2(0:NAR) 
DATA KAY, QFLX, QFLXO,XPO, CONST, Q, N / 0.067,25.0,25.0,0.3,5.4,0,0.667/ 
DATA K, DEN, C, TP, TA, PI / 0.000432, 1190.0, 4.12, 375.0, 25.0, 3.14/ 

DATA TOL, H /0.0001,1 / 


OPEN(UNIT=11, FILE='CASE3.DAT’, STATUS = 'UNKNOWN’) 
SUM1(0) = 0.0 
SUM2(0) = SUM1(0) 
TAU = ((PI*K*DEN*C)*(TP-TA)) / (4*QELXO**2) 
VP(0) = (1/TAU)*((KAY*(Q+(QELX*CONST*XP0))**N)-XPO) 
DO 100 I = 0,999 
M(1) = CONST 
Md+1) = M() 
VP(I+1) = VP(1) 
C..... ITERATION LOOP 


K=0 
A(ly=0.0 
10 K = K+1 

IF(K .GT.NMAX) GOTO 45 
G1 = M(1) * VP() 
G2 = M(I+1) * VP(I+1) 
G3 = VP(I) + VP (I+1) 
SUM1(I+1) = SUM1(1I)+G1+G2 
SUM2(I+1) = SUM2(1)+G3 


VP(I+1) = (1/TAU)*(KAY*(Q+QFLX*(M(Ii+1)*XPO+ 
7 ((H/2)*(SUM1(I+1)+G1+G2))))**(N)-(XPO+((H/2)*(SUM2(I+1)+G3)))) 
X(K) = VP(I+1) 


Are FIRST ITERATION. NO CONVERGENCE TEST 
IF(K .EQ.1) GOTO 10 
ERR = ABS((X(K) - X(K-1))/X(K)) 

IF (ERR .LE. TOL) GOTO 50 

GOTO 10 
45 WRITE(11,200)'DID NOT CONVERGE’, VP(I+1), ERR 
200 FORMAT(AI6, E8.2, 5X, E8.2,) 
50 WRITE(11,201)';CONVERGED AT’, K,'-TH ITERATIONS WITH ERR=', ERR 
201 FORMAT(AI3, 15,A24,3X,E8.2) 

WRITE(11,444) I+1, VP(I+1) 

100 CONTINUE 
444. FORMAT(10X, IS, 5X, E11.6) 

STOP 

END 


Program 4 Exact Solution for n=2/3 
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APPENDIX C VARIABLES AND DATA USED FOR GENERALIZED 
FLAME SPREAD MODEL 
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Material Data 


T.. - initial or ambient temperature 
k - thermal conductivity 
p - density 
c - specific heat 
T._ - ignition Temperature 
ig 


AH - heat of vaporization 
L - heat of gasification 
AH. - heat of combustion 
a - thermal diffusivity 


¢- thickness 


m” - burnable mass per unit area 


igs. 


Properties Name of 
Variable in 


Program 


ae AT 300 


= vial er 
Me 


TABLE C.1 Names of Variables and Data Used for Material 
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Ignitor Characteristic 


Qie - energy output rate 


w - width of wall heated 


Xo ~ ignitor flame height 


At. - duration of burning 
ig 


Q' ig ~ effective energy rate per wall width 


Name of 


Properties 


Variable in 


Program 


: ty i 
1g 


WwW 


AS ol 
en ee i ae 
ee 


TABLE C.2 Names of Variables and Data Used for Ignitor Characteristic 


iff 


Heat Flux 


q" fs heat flux from wall 


q' 6 on heat flux from ignitor 


Properties Name of 


Variable in 


Program 


TABLE C.3 Names of Variables and Data Used for Heat Flux 
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Flame Height 


k, - coefficient 


Nn - power 


Properties Name of 


Variable in 


Program 


0.01, n=1 


0.067, n=2/3 


TABLE C.4 Names of Variables and Data Used for Flame Height 


Te, 


Computational Parameters 


h - time step 


a8 - maximum time 
max 


€ - tolerance for convergence 


o - Stefan Boltzmann constant 


Properties Name of 


Variable in 


Program 


TABLE C.5 Names of Variables and Data Used for Computational Parameters 
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Computed Parameters 


t. - ignition time 
T, —T 
TU co 
t. = — kpc — 
ig 4 q fig 


At,- spread time 


Tl 
Ate = ‘4: kpc 


8. - steady penetration depth 


2kL 
ds = 


ue 4 
c(q fan oT; ) 


& 


q" - net flame heat flux 


4 
g 


ou : 


ey Ay 


* 
t™ -bum time constant 
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Properties Name of 


Variable in 


Program 


TABLE C.6 Names of Variables Used for Computed Parameters 
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Output Variables 


T(i+1) - time after ignition 

13 i) - time of arrival of flame tip at position X=x (i) 

m''(j,1) - burning rate per unit area at time T(j) and position K=X (1) 
T, (a) - burnout time at x=x (i) 

x (a) - burnout position at T=T(1) 


Q(i) - total energy release rate 
x (1) - pyrolysis front position 
X,(i) - flame tip position 


Vp(i) - velocity of the pyrolysis front 
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Properties Name of Units 
Variable in 
Program 
i posal 
: ie —_— 


TABLE C.7 Names of Variables Used for output 
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APPENDIX D THE PROGRAM OF GENERALIZED FLAME SPREAD 
MODEL 
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PARAMETER(NMAX=1005) 
DECLATION PART 
REAL AT 
DATA AT /20.0/ 
MATERIAL DATA 
REAL K1, DEN, C, IGT, HV, L, HC, THDI, THI, BM 
DATA K1, DEN, C, IGT, L /0.000432,1190.0,4.12,180.0,2770.0/ 
DATA HC, THI /25000.0,0.005/ 
IGNITOR CHARACTERISTICS 


REAL QIG, W, XFIG, DELTIG, EQIG 
DATA QIG, W, DELTIG /100.0,0.5,200.0/ 


HEAT FLUX 


REAL QF, QFIG 
DATA QF, QFIG /30.0,30.0/ 


FLAME HEIGHT 


REAL KF, P 
DATA KF, P /0.01,1.0/ 


COMPUTATIONAL PARAMETERS 


REAL H, TAUMAX, SIG 
DATA H, SIG, TAUMAX /1.0,5.67E-11,1000.0/ 


COMPUTED PARAMETERS 

REAL TIG, TF, DS, NETQ, TAUS, DELIG, DEL1 
MAXIMUM NUMBER OF STEPS, N 

REAL N 
FOR ITERATION 


REAL TOL 
DATA TOL /0.05/ 
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Siege 


VARIABLES 


REAL TAU(2000), TAUF(NMAX), TAUP(NMAX), TAUB(NMAX) 
REAL XP(NMAX), XF(NMAX), XB(NMAX) 

REAL M, VP(NMAX) 

REAL Q1, Q2, TQ(NMAX) 

REAL SUM1(0:2000), SUM2(NMAX) 

INTEGER I 


COMMON STATEMENT 


COMMON /TIME/ TAU, TAUP, TAUF 

COMMON /BUNT/ TAUB 

COMMON /HEIGHT/ XB, XP, XF, VP 

COMMON /BRT/ BM, SUM1 

COMMON /DELTA/ DEL 1 

COMMON /RTM/ DS, TOL, NETQ, IGT, AT, HV, TAUS, THDI 


COMMON /SPD/ DELTIG, TIG, XFIG, HC, SUM2, KF, P, TF, TQ, W 


COMMON /CONT/ H, K1 


OPEN(UNIT=11, FILE='TEST3.DAT’, STATUS = 'UNKNOWN') 


INITIALIZATION PART 

MATERIAL DATA 
HV =L- (C*(IGT - AT)) 
THDI = (K1) / (DEN * C) 
BM = DEN*THI 

FLAME HEIGHT 


XFIG = KF * (QIG/W)**P 
EQIG = (XFIG/KF)**(1/P) 


COMPUTED PARAMETERS 
TIG = (3.14159/4) * (K1*DEN*C) * (IGT-AT)/QFIG)**2 
TF = (3.14159) * (K1*DEN*C) * (IGT-AT)/QF)**2 
DS = 2*Ki*L)/ (C * (QF =(SIG*IGI**4))) 
NETQ = QF - (SIG*IGT**4) 
TAUS = (DS**2*HV) / (6*THDI*L) 
MAXIMUM NUMBER OF STEP 


N = TAUMAX/H 
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Cee INITIAL CONDITION 


SUM1(0) = 0.0 
SUM2(1) = 0.0 
=o) 

TAU() = 0.0 
TAUP(I) = 0.0 
TAUF(I) = -TIG 


Cre DEL1 IS DEL IN THE INITIAL CONDITION. 
DELIG = SQRT((6*THDI) * (TAU(D - TAUF())) 
DEL1 = DELIG 
M = (NETQ - (((2*K1)*(IGT-AT)) / (DEL1))) / HV 
Coes FIND BURNOUT TIME (TAUB) 
CALL BURNOUT(I) 
IF((TAUB(1) .GE. TAU(D)) .OR. ((DELTIG-TIG) .GE. TAU(1))) GOTO 10 
EQIG = 0.0 
10 EQIG = EQIG 
Q1=HC*M 
XB(I) = 0.0 
XF(1) = XB(1) + KF * (EQIG + Q1) 
XP(I) = XFIG 
VP(I) = (XF(I) - XP(D) / TF 
Go TOTAL ENERGY RELEASE RATE 
TQ(D = (EQIG + Q1) * W 
WRITE(11,20) TAU(D,',',XB(1),',',. XP(D),',’,XF(),',, VP(D,',", TOM) 


20 FORMAT(F7.2,A2,1X,E11.6, A2,1X,E11.6,A2,1X,E11.6, A2,1X,E11.6, 
A2,1X,E11.6) 
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1 0).0 


30 


MAIN LOOP 

DO 30I1=1,N 
TAU(I+1) = TAU() +H 
IF (XP(1) .LE. XF(1)) THEN 
TAUF(I+1) = 0.0 
ELSE 
TAUF(I+1) = TAUF(D 
END IF 

FIND BURNOUT TIME 

CALL BURNOUT (I+1) 


FIND BURNOUT POSITION 


IF (TAU(I+1) .GE. TAUB(1)) THEN 


CALL SEARCHB(I+1) 
BSE 

XB(I+1) = 0.0 

END IF 


FIND SPREAD RESULT 


CALL SPREAD(I+1) 


FIND TIME FLAME TIP REACHED FIRST 


IF (XP(I+1) .LE. XF(1)) THEN 


TAUF(I+1) =0.0 
ELSE — 
CALL SEARCHF(I+1) 
END IF 


CONTINUE 


STOP. 
END 
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Se 


pee! 


FIND BURNOUT TIME 
SUBROUTINE BURNOUTV(J) 
PARAMETER(NMAX=1005) 


REAL TAU(2000), TAUP(NMAX), TAUF(NMAX) 
REAL TAUB(NMAX) 

REAL BM, SUM1(0:2000) 

REAL DS, TOL, NETQ, IGT, AT, HV, TAUS, THDI 
REAL H, K1 

REAL A, B, G1, G2, G3 

INTEGER I, J, N1 


COMMON /TIME/ TAU, TAUP, TAUF 

COMMON /BUNT/ TAUB 

COMMON /BRT/ BM, SUM1 

COMMON /RTM/ DS, TOL, NETQ, IGT, AT, HV, TAUS, THDI 
COMMON /CONT/ H, K1 


IF(J EQ. 1) THEN 
TAUP(J) = 0.0 
ELSE 

I=J-1 

TAUP(J) = TAU(I+1) 
END IF 


K=0 
K = K+1 
CALL ROOTM(K, J, A) 


A = DEL AT K AND J, B =DEL AT K+1, J 
G1 = M(K,J), G2 = M(K+1,J) 


G1 = (NETQ- (((2*K1)*(IGT-AT)) / A)) / HV 


IF(G1 .LT. 0) THEN 
G1 = 0.0 

ELSE 

G1=Gl 

END IF 


CALL ROOTM(K+1, J, B) 

G2 = (NETQ - (((2*K1)*(IGT-AT)) / B)) / HV 
T(G22 LT? 0) THEN 
G2 = 0.0 


ELSE 
G2 = G2 
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END IF 


SUM1(K) = SUM1(K-1) + G1 + G2 
G3 = (H/2) * SUM1(K) 


IF(G3 .GE. BM) THEN 
Ni=K+1 

TAUB(J) =(N1-1)*H 
ELSE 

GOTO 999 

ENDIE 


RETURN 
END 
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70 
60 
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FIND THERMAL PENETRATION DEPTH 
SUBROUTINE ROOTM(, J, DEL3) 
PARAMETER(NMAX=1005) 


REAL TAU(2000), TAUP(NMAX), TAUF(NMAX) 
REAL DS, TOL, NETQ, IGT, AT, HV, TAUS, THDI 
REAL DEL3, ERR, DELIG, DEL1, X(0:305) 
INTEGER I, J, K 


COMMON /TIME/ TAU, TAUP, TAUF 

COMMON /DELTA/ DEL1 

COMMON /RTM/ DS, TOL, NETQ, IGT, AT, HV, TAUS, THDI 
COMMON /CONT/ H, K1 


IF(I EQ. 1) THEN 
TAU() = 0.0 

ELSE 
TAU() = TAU(-1) +H 
END IF 


DELIG = SQRT((6*THDI) * (TAUP() - TAUF(J))) 
X(0) = DEL1 


ITERATION LOOP TO FIND DELTA 
X(K) IS DELTA 

DO 50 K = 1, 300 

IF(K .GT. 300) GOTO 70 


X(K) =DS-((DS-DELIG)*EXP(((DELIG-X(K-1))/DS)-((TAU(D- 
TAUP(J))/TAUS))) 


IF(X(K) .LT. 0) THEN 
X(K) = X(K-1) 

GOTO 60 

ELSE. 

X(K) = X(K) 

END IF 


ERR = ABS((X(K) - X(K-1)) / X(K)) 


IF(ERR .LE. TOL) GOTO 60 

CONTINUE 
WRITE(11,80) ‘DID NOT CONVERGE’, X(K), ERR 
DEL3 = X(K) 
FORMAT(A16, 2X, E11.6, 5X, E11.6) 
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RETURN 
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Cz FIND J SUCH THAT TAUB(J) = TAU(I) 

SUBROUTINE SEARCHB(I) 

PARAMETER(NMAX=1005) 
REAL TAU(2000), TAUP(NMAX), TAUF(NMAX), TAUB(NMAX) 
REAL XB(NMAX), XP(NMAX), XF(NMAX), VP(NMAX) 
INTEGER I, J 
COMMON /TIME/ TAU, TAUP, TAUF 
COMMON /BUNT/ TAUB 
COMMON /HEIGHT/ XB, XP, XF, VP 

DO 90J=1,1 


IF(TAUB(J) .GE. TAU(I)) GOTO 700 


90 CONTINUE 

700 XB() = XP(J) 
RETURN 
END 
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FIND SPREAD DATA 
SUBROUTINE SPREAD(J) 
PARAMETER(NMAX=1005) 


REAL TAU(2000), TAUP(NMAX), TAUF(NMAX), TAUB(NMAX) 

REAL XB(NMAX), XP(NMAX), XF(NMAX), VP(NMAX) 

REAL DS, TOL, NETQ, IGT, AT, HV, TAUS, THDI 

REAL DELTIG, TIG, XFIG, HC, SUM2(NMAX), SUM3(0:NMAX), KF, P, TE, 
REAL TQ(NMAX), W 

REAL X(0:1005), H, K1, SUM4(0:NMAX), D, E, F, G 

REAL EQIG, ERR1, A, B, C, M1, M2, M3, M4, MS5,FI(1005), FT1(1005), Q1, Q2 
INTEGER I, J, K 


COMMON /TIME/ TAU, TAUP, TAUF 

COMMON /BUNT/ TAUB 

COMMON /HEIGHT/ XB, XP, XF, VP 

COMMON /RTM/ DS, TOL, NETQ, IGT, AT, HV, TAUS, THDI 
COMMON /SPD/ DELTIG, TIG, XFIG, HC, SUM2, KF, P, TF, TQ, W 
COMMON /CONT/ H, K1 


SUM4(0) = 
SUM3(0) = 


FIND EFFECTIVE ENERGY RATE PER WALL WIDTH 


IF(((DELTIG - TIG) .GE. TAU(J)) .OR. (TAUB(1) .GE. TAU(J))) THEN 
EQIG = (XFIG / KF) ** (1/P) 


ELSE 
EQIG = 0.0 
END IF 
FIND Q1 
M1 = M(J,1) 


CALL ROOTMGJ, 1, A) 


IF(TAUB(1) .LT. TAU(J)) THEN 

M1 = 0.0 

ELSE 

M1 = (NETQ - (((2*K1)*(IGT-AT)) / A)) / HV 
END IF 


IF(M1 .LT. 0) THEN 
M1 = 0.0 

ELSE 

M1=M1 

END IF 
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Q1 = HC * M1 * XFIG 
C.... ITERATION LOOP TO FIND VP 

VP(J) = VP(J-1) 

DO 100 K = 1, 1000 


IF (K .GT. 1) GOTO 102 
IF (K .GE. 1000) GOTO 110 
IF (J-2 .EQ. 0) GOTO 102 


C... @BFIND'O? 
DO 101 I= 1,J-2 
C... M2=M(I,J-1), M3 = (MJ, J) 
CALL ROOTM(,, I, B) 


IF(TAUB(]) .LT. TAU(J)) THEN 

M2 = 0.0 

ELSE 

M2 = (NETQ - (((2*K1)*(IGT-AT)) / B)) / HV 
END IF 


IF(M2 .LT. 0) THEN 
M2 = 0.0 

ELSE 

M2 = M2 

END IF 


D = M2 * VP(I) 
CALL ROOTM(, I+1, C) 


IF(TAUB(I+1) .LT. TAU(J)) THEN 

M3 = 0.0 

BESE 

M3 = (NETQ - (((2*K1)*(IGT-AT)) / C)) /HV 
END IF 


IF(M3 .LT. 0) THEN 
M3 = 0.0 

ELSE 

M3 = M3 

END IF 


E = M3 * VP(I+1) 
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SUM3(I) = SUM3(I-1) + VP(I) + VP(I+1) 
SUM4(I) = SUM4(I-1) +D+E 


IF(I .EQ. J-2) GOTO 102 
101 CONTINUE 
102 CALL ROOTMGJ, J-1, F) 


IF(TAUB(-1) .LT. TAU(J)) THEN 

M4 = 0.0 

BES 

M4 = (NETQ - (((2*K1)*(IGT-AT)) / F)) / HV 
END IF 


IF(M4 .LT. 0) THEN 
M4 = 0.0 

ELSE 

M4 = M4 

END IF 


CALL ROOTMGJ, J, G) 


IF(TAUB() .LT. TAU(J)) THEN 

M5 = 0.0 

ELSE 

MS5 = (NETQ - (((2*K1)*(IGT-AT)) / G)) /HV 
END IF 


IF(M5 .LT. 0) THEN 
M5 = 0.0 

ELSE 

M5 = M5 

END IF 


IF(J-2 .EQ. 0) THEN 
SUM3(J-2) = 0.0 
SUM4(J-2) = 0.0 

EUSE 

SUM3(J-2) = SUM3(J-2) 
SUM4(J-2) = SUM4(J-2) 
END IF 


FI(K) = (H/2) * (SUM4(J-2) + (F*VP(J-1)) + (G*VP(J))) 
Q2 = HC * FI(K) 


Ge FIND FLAME HEIGHT(XF) 
XF(J) = ((KF) * (EQIG + Q1] + Q2)**P)) + XBQ) 


Ce. FIND XP 


oa 


Cx 


100 
110 
130 
120 


150 


FI1(K) = (H/2) * (SUM3(J-2) + VP(J)) 
XP(J) = XFIG + FI1(K) 


FIND VP 


VP(J) = (XF(J) - XP(J)) / TF 

X(K) = VP(J) 

ERR1 = ABS((X(K) - X(K-1)) / X(K)) 

IF (ERR1 .LE. TOL) GOTO 120 

CONTINUE 
WRITE(11,130) 'DID NOT CONVERGE’, X(K), ERR1 
FORMAT(AI6, 2X, E11.6, 5X, E11.6) 

TQ(J) = (EQIG + Q1 + Q2)* W 


WRITE(11,150)TAU(J),','",XB(J),',’,XP(J),',',XF(J),',’, VP(),',, TQ) 
FORMAT(F7.2,A2,1X,E11.6,A2,1X,E11.6,A2,1X,E11.6, A2,1X,E11.6, 
A2,1X,E11.6) 


RETURN 
END 
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FIND THE TIME THE FLAME TIP FIRST REACHED X=XP(J) 
SUBROUTINE SEARCHF(J) 
PARAMETER(NMAX=1005) 


REAL TAU(2000), TAUP(NMAX), TAUF(NMAX) 

REAL XB(NMAX), XP(NMAX), XF(NMAX), VP(NMAX) 
REAL H 

INTEGER J, K 


COMMON /TIME/ TAU, TAUP, TAUF 
COMMON /HEIGHT/ XB, XP, XF, VP 
COMMON /CONT/ H 


IF(XF(K) .GT. XP(J)) THEN 
TAUF() = (K-1)*H 
ELSE 

TAUF(J) = 0.0 

GOTO 160 

END IF 


RETURN 
END 
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APPENDIX E THE RESULTS OF GENERALIZED FLAME SPREAD 
MODEL 
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FIGURE E.1_ Flame tip position, pyrolysis front position, and burnout position as a 
function of time of generalized flame spread model for PMMA. 
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FIGURE E.2 Bumout effect of Flame tip position as a function of time of generalized 
flame spread model for PMMA. 
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FIGURE E.3 _ Ignitor effect of Flame tip position as a function of time of generalized 
flame spread model for PMMA. 
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FIGURE E.4 Velocity of the pyrolysis front as a function of time of generalized 
flame spread model for PMMA. 
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FIGURE E.5 Burnout position as a function of time of generalized flame spread model 
for PMMA. 
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FIGURE E.6 Energy release rate as a function of time of generalized flame 
spread model for PMMA. 


106 


APPENDIX F VARIABLES AND DATA USED FOR COMPARISON 
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Properties Name of Data Units 
Variable in 
Program 
2:63"10— kW/m.°C 


1190.0 kg/m3 


2.09 
363.0 


TABLE F.1 Orloff, de Ris, and Markstein‘s Data for Comparison 
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Properties Name of Units 
Variable in 
eee 


= 7a as 


TABLE F.2 LIFT‘s Data for Comparison 
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APPENDIX G PROGRAMS FOR COMPARISON 


REAL KAY, QFLX,QFLX0, XPO, CONST, Q, K, DEN, C, TP, TA 

REAL TAU, XP(0:1000), XF(0:1000), VP(0:1000), M(0:1000) 

DATA KAY, QFLX, QFLXO,XPO, CONST, Q, N / 0.01,25.0,25.0,0.02,9.6,0, 1/ 
DATA K, DEN, C, TP, TA, PI / 0.000263, 1190.0, 2.09, 363.0, 20.0,3.14/ 


OPEN(UNIT=11, FILE='CASE1.DAT', STATUS = 'UNKNOWN') 


TAU = ((PI*K*DEN*C)*(TP-TA)) / (4*QFLXO**2) 

DO 100 I=0,1000 
M() = CONST 
XP(I) = XPO*EXP((KAY*M(1)*QFLX-1)*I/TAU) 
XF(I) = KAY*((Q+(QFLX*M(I)*XP(1)))**N) 
VP(1) = (XF(1) -XP(1))/TAU 

WRITE(11,444) I, XP(1), XF(1), VP(1) 


100 CONTINUE 

444 FORMAT(10X, 15, 5X, E11.6, 5X, E11.6, SX, E11.6) 
STOP 
END 


Program1 Exact Solution for n=1 used Orloff, de Ris, and Markstein‘s Data 
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REAL KAY, QFLX,QFLX0, XPO, CONST, Q, K, DEN, C, TP, TA 
REAL TAU, XP(0:1000), XF(0:1000), VP(0:1000), M(0:1000) 
DATA KAY, QFLX, QFLXO,XPO, CONST, Q, N / 0.067,25.0,25.0,0.02,9.6,0,0.667/ 
DATA K, DEN, C, TP, TA, PI / 0.000263, 1190.0, 2.09, 363.0, 20.0,3.14/ 
OPEN(UNIT=11, FILE='CASE2.DAT’,, STATUS = 'UNKNOWN’) 
TAU = ((PI*K*DEN*C)*(TP-TA)) / (4*QFLXO**2) 
DO 100 I=0,1000 
M(I) = CONST 
A = EXP((-1*1)/(3*TAU)) 
B = 1 - (XPO**(1,/3.) (KAY *((M(1)*QFLX)**(2./3.))) 
C=A*B 
D=(1-C)**3 
E = (KAY**3)*(M(1)*QFLX)**2 
XP(1) = D*E 
XF(I) = KAY*((Q+(QFLX*M(I)*XP(D))**N) 
VP(I) = (XF() -XP(D)/TAU 
WRITE(11,444) I, XP), XF(I), VP(1) 


100 CONTINUE 

444 FORMAT(1X,15,1X,E11.6,1X,E11.6,1X,E11.6) 
STOP 
END 


Program 2 _ Exact Solution for n=2/3 used Orloff, de Ris, and Markstein‘s Data 


Abe 


REAL KAY, QFLX,QFLX0, XPO, CONST, Q, K, DEN, C, TP, TA 

REAL TAU, XP(0:1000), XF(0:1000), VP(0:1000), M(0:1000) 

DATA KAY, QFLX, QFLXO,XPO, CONST, Q, N / 0.01,25.0,25.0,0.02,9.6,0, 1/ 
DATA K, DEN, C, TP, TA, PI / 0.000346, 1180.0, 2.5, 363.0, 20.0,3.14/ 


OPEN(UNIT=11, FILE='CASE1.DAT’, STATUS = 'UNKNOWN') 


TAU = ((PI*K*DEN*C)*(TP-TA)) / (4*QFLXO**2) 

DO 100 I=0,1000 
M(I) = CONST 
XP(I) = XPO*EXP((KAY*M(1)*QFLX-1)*I/TAU) 
XF(I) = KAY*((Q+(QFLX*M(I)*XP(1))**N) 
VP(I) = (XF() -XP()/TAU 

WRITE(1 1,444) I, XP(I), XF(1), VP(1) 


100 CONTINUE 

444 FORMAT(10X, I5, 5X, E11.6, 5X, E11.6, 5X, E11.6) 
STOP 
END 


Program 3 Exact Solution for n=1 used LIFT Data 
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REAL ‘KAY, ‘QFLX,OFUX0, XPO, ‘CONST, QO) Ko DEN, COEPALTA 
REAL TAU, XP(0:1000), XF(0:1000), VP(0:1000), M(0:1000) 


DATA KAY, QFLX, QFLXO,XPO, CONST, Q, N / 0.067,25.0,25.0,0.02,9.6,0,0.667/ 
DATA K, DEN, C, TP, TA, PI / 0.000346, 1180.0, 2.5, 363.0, 20.0,3.14/ 


OPEN(UNIT=11, FILE='CASE2.DAT’, STATUS = 'UNKNOWN’) 
TAU = ((PI*K*DEN*C)*(TP-TA)) / (4*QFLXO**2) 
DO 100 I=0,1000 
M(1) = CONST 
A = EXP((-1*1/(3*TAU)) 
B = 1 - (XPO**(1./3.))(KAY*((M()*QFLX)**(2./3.))) 
C=A*B 
D=(1 -C)**3 
E = (KAY**3)*(M(1)*QFLX)**2 
XP(I) = D*E 
XF(I) = KAY*((Q+(QFLX*M(1)*XP(1)))**N) 
VP(I) = (XF(I) -XP(1))/TAU 
WRITE(11,444) I, XP(1), XF(1), VP(1) 


100 CONTINUE 

444 FORMAT(1X,15,1X,E11.6,1X,E11.6,1X,E11.6) 
STOP 
END 


Program4 Exact Solution for n=2/3 used LIFT Data 
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